Constrained optimization approach to compander design for OFDM PAPR reduction

ABSTRACT

In the method for constrained optimization for compander design for OFDM PAPR, wherein the improvement comprises the step of converting a Rayleigh amplitude distribution from which superior performing companders are derived, whereby a constrained optimization problem may be solved using Lagrange multipliers.

RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Application No. 62/036,388, filed Aug. 12, 2014. This application is herein incorporated by reference in its entirety for all purposes.

FIELD OF THE INVENTION

The invention relates to signal processing and, more particularly, to companders for reducing peak to average power ratio (PAPR) in orthogonal frequency division multiplexing (OFDM) signals.

BACKGROUND OF THE INVENTION

Orthogonal frequency-division multiplexing (OFDM) is a method of encoding digital data on multiple carrier frequencies. OFDM has developed into a popular scheme for wideband digital communication, used in applications such as digital television and audio broadcasting, DSL Internet access, wireless networks, powerline networks, and 4G mobile communications.

For a sufficiently large number of tones, OFDM signal values are Gaussian, or normally, distributed, by the central limit theorem, which states that, given certain conditions, the arithmetic mean of a sufficiently large number of iterates of independent random variables, each with a well-defined expected value and well-defined variance, will be approximately normally distributed, regardless of the underlying distribution.

The Gaussian distributed signal components generate Rayleigh distributed signal amplitude values, Rayleigh distributions being continuous probability distributions for positive-valued random variables, specifically the distribution of the magnitude of a two-dimensional random vector whose coordinates are independent, identically distributed, mean 0 normal variables. In signals typical of the current state of the art, the Rayleigh distribution is typically heavy-tailed, i.e. a larger proportion of its population rests within its tail than would under a normal distribution. The long tail of the Rayleigh distribution, i.e. the portion of the distribution having a large number of occurrences far from the “head” or central part of the distribution, indicates the presence of large amplitude values, resulting from constructive interference. These large amplitude values can saturate the power amplifier, used in the processing chain, thus creating signal distortion which reduces demodulation performance.

This causes existing techniques for Orthogonal Frequency-Division Multiplexing (OFDM) to suffer from large peak-to-average power ratio (PAPR). The PAPR is the peak amplitude squared (giving the peak power) divided by the mean square value squared (giving the average power).

Many techniques have been developed to mitigate the PAPR problem, including signal companding. Companders are attractive because of their simplicity and effectiveness. Companders apply an amplitude weight to the signal, generally downweighting large amplitude values, while upweighting smaller values to maintain constant power. The large amplitude downweighting, while simultaneously maintaining the same power level, reduces the PAPR. The amplitude weighting introduces signal distortion, so it is important to design companders that introduce minimal distortion to limit demodulation errors, while also providing out-of-band power rejection.

Early solutions adapted companders from speech and audio coding, such as μ-law and A-law companders. Subsequent approaches decomposed the amplitude interval into disjoint regions over which different amplitude weighting functions were used. Modern companders are derived by transforming the Rayleigh amplitude distribution into a more favorable distribution to reduce PAPR. Among these transformations, piecewise transformations, i.e. a function composed of straight-line sections, have been used, either uniform, linear, nonlinear, or piecewise linear.

A need still exists, however, for solutions which minimally alter the original Rayleigh amplitude distribution.

SUMMARY OF THE INVENTION

One embodiment of the present invention provides a signal transmitting apparatus for an orthogonal frequency division multiplexing (OFDM) system, comprising: a compander configured to compress and expand a transmitted signal wherein the compander comprises a Rayleigh amplitude distribution to a constrained optimization problem which may be solved using Lagrange multipliers.

One embodiment of the present invention provides a signal transmitting apparatus for orthogonal frequency division multiplexing (OFDM) system, comprising: a companding feedback module configured to compress and expand a transmitted signal; wherein the companding feedback module comprises a Rayleigh amplitude distribution to a constrained optimization problem which may be solved using Lagrange multipliers.

One embodiment of the present invention provides a method of deriving a compander for orthogonal frequency division multiplexing comprising: converting a Rayleigh amplitude distribution to a constrained optimization problem which may be solved using Lagrange multipliers.

Another embodiment of the present invention provides such a method wherein the Rayleigh amplitude distribution is decomposed into disjoint regions over which a plurality of amplitude weighting functions are used.

A further embodiment of the present invention provides such a method wherein the constrained optimization approach comprises optimizing an objective function with respect to changeable elements in the presence of constraints on those elements.

Yet another embodiment of the present invention provides such a method wherein the Rayleigh amplitude distribution is a probability distribution function and a first the constraint is that the Rayleigh amplitude distribution must integrate to unity.

A yet further embodiment of the present invention provides such a method wherein a second the constraint is that power must be constant across the distribution.

Still another embodiment of the present invention provides such a method wherein an optimal solution is minimally perturbed to generate only non-negative solutions.

A still further embodiment of the present invention provides such a method of claim 5 wherein the Rayleigh amplitude distribution is defined as:

${f_{R}(x)} = {\frac{2x}{\sigma^{2}}{{\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2\;}}}.}}$

Even another embodiment of the present invention provides such a method wherein further comprising determining a function g(x) which minimizes:

∫₀^(A)(f_(R)(x) − g(x))²𝕕x.

An even further embodiment of the present invention provides such a method wherein g(x) is chosen to be of piecewise linear form.

A still even another embodiment of the present invention provides such a method wherein g(x) is defined as

${g(x)} = {\sum\limits_{i = 1}^{N}{{U_{i}(x)}.}}$

A still even further embodiment of the present invention provides such a method wherein the dynamics of the system are defined as: Λ(β₁, . . . β_(N+1), λ₁, λ₂)=F(β₁, . . . β_(N+1))+λ₁[g₁(β₁, . . . β_(N+1))−1]+λ₂[g₂(β₁, . . . β_(N+1))−σ²].

Still yet another embodiment of the present invention provides such a method wherein a first constraint, that the cumulative distribution function must integrate to unity, is expressed as: g₁(β₁, . . . β_(N+1))=1.

A still yet further embodiment of the present invention provides such a method wherein a second constraint, that power must be constant, is expressed as: g₂(β₁, . . . β_(N+1))=σ².

Even yet another embodiment of the present invention provides such a method wherein the constrained optimization problem can be defined as: Minimize:

F(β₁, …  β_(N + 1)) = ∫₀^(A)(f_(R)(x) − g(x))²𝕕x for β_(i), subject to

∫₀^(A)g(x)𝕕x = 1  and  ∫₀^(A)x²g(x)𝕕x = ∫₀^(∞)x²f_(R)(x)𝕕x = σ², using Lagrange multipliers.

An even yet further embodiment of the present invention provides such a method wherein the compander is defined as:

${F_{y}^{- 1}{F_{R}(x)}} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$ for a_(i)≦x<a_(i+1), wherein

$Z_{i} = {Z_{i - 1} + {\frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}.}}$

Still even yet another embodiment of the present invention provides such a method wherein a decompander for use with the compander is defined as:

${x = {{{sgn}(x)}\sqrt{{- \sigma^{2}}{\ln\left( {{- {\gamma\left( {y + \frac{\delta}{2\gamma}} \right)}^{2}} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}} \right)}}}},$ where the sgn function refers to retaining either the positive or negative square-root, depending on the sign of x; and wherein the range of the compander, when a_(i)≦x<a_(i+1), then y_(MIN)≦y<y_(MAX) where: y_(MIN)=min{y(a_(i)),y(a_(i+1))}, y_(MAX)=max{y(a_(i)),y(a_(i+1))}, can be defined as

${y\left( a_{i + 1} \right)} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{a_{i + 1}^{2}}{\sigma^{2\;}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - {\frac{\delta}{2\gamma}.}}$

A still even yet further embodiment of the present invention provides such a method wherein a Hessian Matrix is used to determine critical point types.

The features and advantages described herein are not all-inclusive and, in particular, many additional features and advantages will be apparent to one of ordinary skill in the art in view of the drawings, specification, and claims. Moreover, it should be noted that the language used in the specification has been principally selected for readability and instructional purposes, and not to limit the scope of the inventive subject matter.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a graph showing piecewise linear approximation to a preferred embodiment of the invention;

FIG. 2 is a graph showing optimal solution, A=1.1;

FIG. 3 is a graph showing optimal solution, A=1.2;

FIG. 4 is a graph showing optimal solution, A=1.3;

FIG. 5 is a graph showing optimal solution, A=1.4;

FIG. 6 is a graph showing optimal solution, A=1.5;

FIG. 7 is a graph showing optimal solution, A=1.6;

FIG. 8 is a graph showing optimal solution, A=1.7;

FIG. 9 is a graph showing optimal solution, A=1.8;

FIG. 10 is a graph showing positive perturbed solution; A=1.4;

FIG. 11 is a graph showing positive perturbed function;

FIG. 12 is a graph showing perturbation solution; A=1.4, k=5;

FIG. 13 is a graph showing perturbation solution; A=1.4, k=6;

FIG. 14 is a graph showing perturbation solution; A=1.4, k=7;

FIG. 15 is a graph showing perturbation solution; A=1.4, k=8;

FIG. 16 is a graph showing symbol error rate performance;

FIG. 17 is a graph showing PAPR reduction performance;

FIG. 18 is a graph showing out-of-band power rejection; and

FIG. 19 is a graph showing sample Lagrange Companders.

DETAILED DESCRIPTION

The present disclosure describes new companders that use a constrained optimization approach for reducing peak-to-average power ratio (PAPR) in orthogonal frequency-division multiplexing (OFDM) signals. A constrained optimization approach, in this context, involves optimizing an objective function, specifically an energy function, with respect to changeable elements, i.e. variables, in the presence of constraints on those variables.

Companders, devices which use both compression and expansion to improve dynamic range and signal-to-noise ratio, in accordance with the present disclosure provide symbol error rate performance improvements over current state-of-the art companders. In digital communications, symbol rate (also known as baud or modulation rate) is the number of symbol changes (waveform changes or signalling events) made to the transmission medium per second using a digitally modulated signal. The symbol rate is measured in baud (Bd) or symbols per second. Each symbol can represent or convey one or several bits of data. The symbol error rate is simply the rate at which signalling events fail to convey the intended data.

The newly-designed companders also provide design flexibility, thereby expanding the space of tradeoffs between demodulation performance, PAPR reduction, and out-of-band power rejection. Furthermore, the new companders provide solutions in operating regions where certain current companders fail to exit; solutions may be derived for cutoff amplitude values that are unobtainable using other companders. The new compander solutions enable tuning the cutoff amplitude value based on power amplifier bandwidth.

Below, we formulate the constrained optimization problem and derive the compander and decompander forms. Through numerical simulation, we generate performance results demonstrating the capability of the new companders.

1. COMPANDER DESIGN USING CONSTRAINED OPTIMIZATION

This section describes the constrained optimization problem and solution, and derives the compander and decompander.

1.1. Constrained Optimization Problem Definition

Denote the Rayleigh amplitude distribution by ƒ_(R)(x) for:

${f_{R}(x)} = {\frac{2x}{\sigma^{2}}{\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}}}$

In the above equation, ƒ_(R)(x) represents the probability density function of the amplitude values of the OFDM signal. This is a statistical description of the spread of amplitude values of the OFDM signal, which follow a Rayleigh distribution. We then look for a function g(x) which minimizes:

$\begin{matrix} {\int_{0}^{A}{\left( {{f_{R}(x)} - {g(x)}} \right)^{2}{\mathbb{d}x}}} & (1) \end{matrix}$

This function is used because we are making the assumption, or operating under the hypothesis, that minimal distortion of the natural signal distribution (i.e., the Rayleigh distribution) will incur minimal signal distortion of the companded function and hence reduce the demodulation errors (i.e., reduce the symbol error rate).

We then choose g(x) to be of piecewise linear form:

${g(x)} = {\sum\limits_{i = 1}^{N}{U_{i}(x)}}$

In the above equation, the piecewise linear form was chosen since it extends previous work on using single linear components, with the view that the increased flexibility of the piecewise approach will lead to better performing companders. Another reason it was used is that it allows a closed form expression for the compander weights to be found, by solving a linear system of equations.

Still referring to the above equation, each U_(i)(x) is a compactly-supported linear segment, in this context, meaning a line segment of finite extent, i.e., of finite length. The U_(i)(x) is defined by equation (2) and for each index i, represents the equation for one of the linear segments shown in FIG. 1.

$\begin{matrix} {{U_{i}(x)} = \left\{ \begin{matrix} {{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)\left( {x - a_{i\;}} \right)} + \beta_{i}} & {{{for}\mspace{14mu} x} \in \left\lbrack {a_{i},a_{i + 1}} \right\rbrack} \\ 0 & {otherwise} \end{matrix} \right.} & (2) \end{matrix}$

and (a₁, . . . a_(N+1)) is a partition of the interval [0, A] and the β_(i) are the ordinates at x=a_(i) (FIG. 1). This interval is a standard interval used in PAPR reduction techniques, and represents an interval of amplitude values, covering amplitude values between 0 and A. Any amplitude value in the original OFDM signal that is larger than A gets mapped to something less than or equal to A, so A is the maximum value the companded value can have.

With the assumed form of g(x), we want to minimize:

$\begin{matrix} {{F\left( {\beta_{1},{\ldots\mspace{14mu}\beta_{N + 1}}} \right)} = {\int_{0}^{A}{\left( {{f_{R}(x)} - {g(x)}} \right)^{2}\ {\mathbb{d}x}}}} & (3) \end{matrix}$

as a function of the β_(i). The function that we are minimizing represents a measure of the deviation of the new, companded function amplitude value distribution from the original OFDM Rayleigh amplitude distribution. The beta values are the values of the piecewise linear component functions at the endpoints, as shown in FIG. 1. (i.e., the ordinates), and described by the functions U_(i)(x).

There are two constraints on g(x); g(x) is a probability density function and so must integrate to unity value over xε[0, A]. This is because, by definition, a probability density function gives an amount of probability per unit of something, in this case, per unit of signal amplitude value. Since the total probability of all of the possibilities adds up (integrates up) to one, then the integral of the probability density function over all of its values must equal the total probability value of one.

$\begin{matrix} {{\int_{0}^{A}{{g(x)}\ {\mathbb{d}x}}} = 1} & (4) \end{matrix}$

This is the unity cumulative distribution function (c.d.f.) constraint. The c.d.f. itself describes the probability that a real-valued random variable x with a given probability distribution will be found to have a value less than or equal to x. In the case of a continuous distribution, it gives the area under the probability density function (p.d.f.) from minus infinity to x.

The second constraint is the constant power constraint:

$\begin{matrix} {{\int_{0}^{A}{x^{2}{g(x)}\ {\mathbb{d}x}}} = {{\int_{0}^{\infty}{x^{2}{f_{R}(x)}\ {\mathbb{d}x}}} = \sigma^{2}}} & (5) \end{matrix}$

In the above equation, equation 5, the integral on the right represents the total power across all possible amplitude values, and its value is equal to the constant value sigma times sigma, hence the total power is the constant value sigma times sigma.

There is actually a third constraint on g(x); because g(x) is a p.d.f., we have a non-negativity constraint: g(x)≧0  (6)

However, we forgo the non-negativity constraint, because, as is discussed later, we can slightly perturb the optimal solution to generate non-negative solutions. The constrained optimization problem then becomes: minimize (3) for the β_(i) subject to (4) and (5). We solve this using Lagrange multipliers, a strategy for finding the local maxima and minima of a function subject to equality constraints.

1.2. Calculation of Partial Derivatives

To solve the constrained optimization problem, we need to calculate the partial derivatives of F. First consider:

$\begin{matrix} {{F\left( {\beta_{1},{\ldots\mspace{14mu}\beta_{N + 1}}} \right)} = {\int_{0}^{A}{\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)^{2}\ {\mathbb{d}x}}}} & (7) \end{matrix}$

This equation results from eq. (1), where we use the piecewise linear form for function g, and call the whole thing F, and in F we explicitly identify the fact that the expression is now a function of the beta parameters, as shown from eq 2.

Because the integrand in (7) is continuously differentiable and integrable for x≠β_(i), then differentiation is permitted under the integral sign and:

$\begin{matrix} {\frac{\partial F}{\partial\beta_{j}} = {{\int_{0}^{A}{\frac{\partial}{\partial\beta_{j}}\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)^{2}\ {\mathbb{d}x}}} = {{\int_{0}^{A}{2\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)\left( {{- \ \frac{\partial}{\partial\beta_{j}}}{\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right){\mathbb{d}x}}} = {{\int_{0}^{A}{2\left( {{f_{R}(x)} - {\sum\limits_{i = 1}^{N}{U_{i}(x)}}} \right)\left( {- \ {\sum\limits_{i = 1}^{N}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}}} \right){\mathbb{d}x}}} = {{{2{\int_{0}^{A}{{- {f_{R}(x)}}{\sum\limits_{i = 1}^{N}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}}}}} + {\sum\limits_{k = 1}^{N}{{U_{k}(x)}{\sum\limits_{i = 1}^{N}{\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}}} = {{{{- 2}{\int_{0}^{A}{{f_{R}(x)}{\sum\limits_{i = 1}^{N}{\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}}} + {2{\int_{0}^{A}{\sum\limits_{k = 1}^{N}{{U_{k}(x)}{\sum\limits_{i = 1}^{N}{\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}}}}} = {{{{- 2}{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{f_{R}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}} + {2{\sum\limits_{i = 1}^{N}{\int_{0}^{A}{{U_{i}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}}} = {{{- 2}{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{f_{R}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}} + {2{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{U_{i}(x)}\frac{\partial{U_{i}(x)}}{\partial\beta_{j}}{\mathbb{d}x}}}}}}}}}}}}} & (8) \end{matrix}$

Now from (2), U_(i)(x) can be rewritten as:

${U_{i}(x)} = {{\left( \frac{x - a_{i}}{a_{i + 1} - a_{i}} \right)\beta_{i + 1}} + {\left( {1 - \frac{x - a_{i}}{a_{i + 1} - a_{i}}} \right)\beta_{i}}}$

so that:

$\begin{matrix} {\frac{\partial U_{i}}{\partial\beta_{j}} = {{\left( \frac{x - a_{i}}{a_{i + 1} - a_{i}} \right)\frac{\partial\beta_{i + 1}}{\partial\beta_{j}}} + {\left( {1 - \frac{x - a_{i}}{a_{i + 1} - a_{i}}} \right)\frac{\partial\beta_{i}}{\partial\beta_{j}}}}} & (9) \end{matrix}$

In the above equation, j represents an index on the beta parameters. For j=1, equation (9), once simplified, becomes:

$\begin{matrix} {\frac{\partial U_{i}}{\partial\beta_{j}} = \left\{ \begin{matrix} {1 - \frac{x - a_{1}}{a_{2} - a_{1}}} & {{{if}\mspace{14mu} i} = 1} \\ 0 & {otherwise} \end{matrix} \right.} & (10) \end{matrix}$

For 1<j≦N, (9) becomes

$\begin{matrix} {\frac{\partial U_{i}}{\partial\beta_{j}} = \left\{ \begin{matrix} {1 - \frac{x - a_{j}}{a_{j + 1} - a_{j}}} & {{{if}\mspace{14mu} i} = j} \\ \frac{x - a_{j - 1}}{a_{j} - a_{j - 1}} & {{{if}\mspace{14mu} i} = {j - 1}} \\ 0 & {otherwise} \end{matrix} \right.} & (11) \end{matrix}$

For j=N+1, (9) becomes

$\begin{matrix} {\frac{\partial U_{i}}{\partial\beta_{j}} = \left\{ {\begin{matrix} \frac{x - a_{N}}{a_{N + 1} - a_{N}} & {{{if}\mspace{14mu} i} = N} \\ 0 & {otherwise} \end{matrix}.} \right.} & (12) \end{matrix}$

So, from (8), and (10), for j=1

$\begin{matrix} {\frac{\partial F}{\partial\beta_{1}} = {{{- 2}{\int_{a_{1}}^{a_{2}}{{f_{R}(x)}\left( {1 - \frac{x - a_{1}}{a_{2} - a_{1}}} \right)\ {\mathbb{d}x}}}} + {2{\int_{a_{1}}^{a_{2}}{{U_{1}(x)}\ \left( {1 - \frac{x - a_{1}}{a_{2} - a_{1}}} \right){\mathbb{d}x}}}}}} & (13) \end{matrix}$

From (8), and (11), for 1<j≦N

$\begin{matrix} {\frac{\partial F}{\partial\beta_{j}} = {{{- 2}{\int_{a_{j - 1}}^{a_{j}}{{f_{R}(x)}\left( \frac{x - a_{j - 1}}{a_{j} - a_{j - 1}} \right)\ {\mathbb{d}x}}}} - {2{\int_{a_{j}}^{a_{j + 1}}{{f_{R}(x)}\ \left( {1 - \frac{x - a_{j}}{a_{j + 1} - a_{j}}} \right){\mathbb{d}x}}}} + {2{\int_{a_{j - 1}}^{a_{j}}{{U_{j - 1}(x)}\left( \frac{x - a_{j - 1}}{a_{j} - a_{j - 1}} \right)\ {\mathbb{d}x}}}} + {2{\int_{a_{j}}^{a_{j + 1}}{{U_{j}(x)}\ \left( {1 - \frac{x - a_{j}}{a_{j + 1} - a_{j}}} \right){\mathbb{d}x}}}}}} & (14) \end{matrix}$

From (8), and (12), for j=N+1

$\begin{matrix} {\frac{\partial F}{\partial\beta_{N + 1}} = {{{- 2}{\int_{a_{N}}^{a_{N + 1}}{{f_{R}(x)}\left( \frac{x - a_{N}}{a_{N + 1} - a_{N}} \right)\ {\mathbb{d}x}}}} + {2{\int_{a_{N}}^{a_{N + 1}}{{U_{N}(x)}\ \left( \frac{x - a_{N}}{a_{N + 1} - a_{N}} \right){\mathbb{d}x}}}}}} & (15) \end{matrix}$

Inspection of (13)-(15) show that the following integral forms are needed:

${K_{R}^{1}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{f_{R}(x)}{\mathbb{d}x}}}$ ${K_{R}^{2}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{{xf}_{R}(x)}{\mathbb{d}x}}}$ ${K_{U}^{1}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{U_{j}(x)}{\mathbb{d}x}}}$ ${K_{U}^{2}\left( {c,d} \right)} = {\int\limits_{c}^{d}{{{xU}_{j}(x)}{\mathbb{d}x}}}$

In the above equations, K is used as a compact way to represent the integral on the right side of the equation. After much simplification, straightforward integration produces:

$\begin{matrix} {\mspace{79mu}{{\int\limits_{c}^{d}{{f_{R}(x)}{\mathbb{d}x}}} = {{\mathbb{e}}^{- \frac{c^{2}}{\sigma^{2}}} - {\mathbb{e}}^{- \frac{d^{2}}{\sigma^{2}}}}}} & (16) \\ {\mspace{79mu}{{\int\limits_{c}^{d}{{{xf}_{R}(x)}{\mathbb{d}x}}} = {{c\;{\mathbb{e}}^{- \frac{c^{2}}{\sigma^{2}}}} - {d\;{\mathbb{e}}^{- \frac{d^{2}}{\sigma^{2}}}} + {\frac{\sigma\sqrt{\pi}}{2}\left\lbrack {{\Phi\left( \frac{d}{\sigma} \right)} - {\Phi\left( \frac{c}{\sigma} \right)}} \right\rbrack}}}} & (17) \\ {{\int\limits_{c}^{d}{{U_{j}(x)}{\mathbb{d}x}}} = {{\left\lbrack {\frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} + \left( {d - c} \right)} \right\rbrack\beta_{j}} + {\left\lbrack {\frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}}} \right\rbrack\beta_{j + 1}}}} & (18) \\ {{\int\limits_{c}^{d}{{{xU}_{j}(x)}{\mathbb{d}x}}} = {{\left\lbrack {\frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} + \left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)} \right\rbrack\beta_{j}} + {\left\lbrack {\frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}}} \right\rbrack\beta_{j + 1}}}} & (19) \end{matrix}$

where Φ(x) denotes the Gaussian error function, a special function (non-elementary) of sigmoid shape that occurs in probability, statistics, and other disciplines:

${\Phi(x)} = {\int\limits_{0}^{x}{\frac{2}{\sqrt{\pi}}{\mathbb{e}}^{- \xi^{2}}{\mathbb{d}\xi}}}$

From (18) and (19) we have:

$\begin{matrix} {K_{U}^{1} = {{{\left\lbrack {\frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} + \left( {d - c} \right)} \right\rbrack\beta_{j}} + {\left\lbrack {\frac{\frac{d^{2}}{2} - \frac{c^{2}}{2}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {d - c} \right)}{a_{j + 1} - a_{j}}} \right\rbrack\beta_{j + 1}}} = {{{M_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}\beta_{j}} + {{N_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}\beta_{j + 1}}}}} & (20) \\ {K_{U}^{2} = {{{\left\lbrack {\frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}} - \frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} + \left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)} \right\rbrack\beta_{j}} + {\left\lbrack {\frac{\frac{d^{3}}{3} - \frac{c^{3}}{3}}{a_{j + 1} - a_{j}} - \frac{a_{j}\left( {\frac{d^{2}}{2} - \frac{c^{2}}{2}} \right)}{a_{j + 1} - a_{j}}} \right\rbrack\beta_{j + 1}}} = {{{M_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}\beta_{j}} + {{N_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}\beta_{j + 1}}}}} & (21) \end{matrix}$

Using (16), (17), (20), and (21) in (13) through (15) after much simplification produces for j=1:

$\begin{matrix} {\frac{\partial F}{\partial\beta_{1}} = {{P_{1}\beta_{1}} + {P_{2}\beta_{2}} + P_{0}}} & (22) \end{matrix}$

Where:

$P_{1} = {{2{M_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2}{a_{2} - a_{1}} \right){M_{1}^{U\; 2}\left( {a_{1},a_{2}} \right)}} + {\left( \frac{2a_{1}}{a_{2} - a_{1}} \right){M_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}}}$ $P_{2} = {{2{N_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2}{a_{2} - a_{1}} \right){N_{1}^{U\; 2}\left( {a_{1},a_{2}} \right)}} + {\left( \frac{2a_{1}}{a_{2} - a_{1}} \right){N_{1}^{U\; 1}\left( {a_{1},a_{2}} \right)}}}$ $\mspace{79mu}{P_{0} = {{{- 2}{K_{R}^{1}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2}{a_{2} - a_{1}} \right){K_{R}^{2}\left( {a_{1},a_{2}} \right)}} - {\left( \frac{2a_{1}}{a_{2} - a_{1}} \right){{K_{R}^{1}\left( {a_{1},a_{2}} \right)}.}}}}$

In equation 22, P is used as a convenient and compact way to represent the expressions shown above.

For 1<j≦N:

$\frac{\partial F}{\partial\beta_{j}} = {{Q_{0}(j)} + {{Q_{1}(j)}\beta_{j - 1}} + {{Q_{2}(j)}\beta_{j}} + {{Q_{3}(j)}\beta_{j + 1}}}$

Where:

${Q_{0}(j)} = {{\left( \frac{- 2}{a_{j} - a_{j - 1}} \right){K_{R}^{2}\left( {a_{j - 1},a_{j}} \right)}} + {\left( \frac{2a_{j - 1}}{a_{j} - a_{j - 1}} \right){K_{R}^{1}\left( {a_{j - 1},a_{j}} \right)}} - {2{K_{R}^{1}\left( {a_{j},a_{j + 1}} \right)}} + {\left( \frac{2}{a_{j + 1} - a_{j}} \right){K_{R}^{2}\left( {a_{j},a_{j + 1}} \right)}} - {\left( \frac{2a_{j}}{a_{j + 1} - a_{j}} \right){K_{R}^{1}\left( {a_{j},a_{j + 1}} \right)}}}$ $\mspace{79mu}{{Q_{1}(j)} = {\frac{2{M_{j - 1}^{U\; 2}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}} - \frac{2a_{j - 1}{M_{j - 1}^{U\; 1}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}}}}$ ${Q_{2}(j)} = {\frac{2{N_{j - 1}^{U\; 2}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}} - \frac{2a_{j - 1}{N_{j - 1}^{U\; 1}\left( {a_{j - 1},a_{j}} \right)}}{a_{j} - a_{j - 1}} + {2{M_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}} - \frac{2{M_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}} + \frac{2a_{j}{M_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}}}$ $\mspace{79mu}{{Q_{3}(j)} = {{2{N_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}} - \frac{2{N_{j}^{U\; 2}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}} + \frac{2a_{j}{N_{j}^{U\; 1}\left( {a_{j},a_{j + 1}} \right)}}{a_{j + 1} - a_{j}}}}$

In the preceding equation for 1<j≦N, Q is used as a convenient and compact way to represent the expressions shown above.

M, as used above, is similarly a convenient, compact expression, seen in eq. 20. It is the term inside the first set of brackets in eq. 20.

N, as used above, is also a convenient, compact expression, seen in eq. 20. It is the term inside the second set of brackets in eq. 20.

For j=N+1:

$\frac{\partial F}{\partial\beta_{N + 1}} = {{R_{1}\beta_{N}} + {R_{2}\beta_{N + 1}} + R_{0}}$

Where:

$R_{0} = {\frac{{- 2}{K_{R}^{2}\left( {a_{N},a_{N + 1}} \right)}}{a_{N + 1} - a_{N}} + \frac{2a_{N}{K_{R}^{1}\left( {a_{N},a_{N + 1}} \right)}}{a_{N + 1} - a_{N}}}$ $R_{1} = {{\left( \frac{2}{a_{N + 1} - a_{N}} \right){M_{N}^{U\; 2}\left( {a_{N},a_{N + 1}} \right)}} - {\left( \frac{2a_{N}}{a_{N + 1} - a_{N}} \right){M_{N}^{U\; 1}\left( {a_{N},a_{N + 1}} \right)}}}$ $R_{2} = {{\left( \frac{2}{a_{N + 1} - a_{N}} \right){N_{N}^{U\; 2}\left( {a_{N},a_{N + 1}} \right)}} - {\left( \frac{2a_{N}}{a_{N + 1} - a_{N}} \right){N_{N}^{U\; 1}\left( {a_{N},a_{N + 1}} \right)}}}$

The R variables used in the preceding paragraphs are just convenient, compact expressions for the quantities shown directly above this paragraph.

To summarize:

$\begin{matrix} {\frac{\partial F}{\partial\beta_{j}} = \left\{ \begin{matrix} {{P_{1}\beta_{1}} + {P_{2}\beta_{2}} + P_{0}} & {{{for}\mspace{14mu} j} = 1} \\ {{Q_{0}(j)} + {{Q_{1}(j)}\beta_{j - 1}} + {{Q_{2}(j)}\beta_{j}} + {{Q_{3}(j)}\beta_{j + 1}}} & {{{for}\mspace{14mu} 1} < j \leq N} \\ {{R_{1}\beta_{N}} + {R_{2}\beta_{N + 1}} + R_{0}} & {{{for}\mspace{14mu} j} = {N + 1}} \end{matrix} \right.} & (23) \end{matrix}$ 1.3. Unity C.D.F. Constraint

From (4), we have:

$\begin{matrix} {{\int\limits_{0}^{A}{{g(x)}{\mathbb{d}x}}} = {{\int\limits_{0}^{A}{\sum\limits_{i = 1}^{N}{{U_{i}(x)}{\mathbb{d}x}}}} = {{\sum\limits_{i = 1}^{N}{\int\limits_{0}^{A}{{U_{i}(x)}{\mathbb{d}x}}}} = {{{\sum\limits_{i = 1}^{N}{\int\limits_{a_{i}}^{a_{i + 1}}{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)\left( {x - a_{i}} \right)}}} + {\beta_{i}{\mathbb{d}x}}} = {{\sum\limits_{i = 1}^{N}{\int\limits_{a_{i}}^{a_{i + 1}}{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)x}}} + {\left\lbrack {\beta_{i} - \frac{a_{i}\left( {\beta_{i + 1} - \beta_{i}} \right)}{a_{i + 1} - a_{i}}} \right\rbrack{\mathbb{d}x}}}}}}} & (24) \end{matrix}$

Let:

${m_{i} = \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}}},{b_{i} = {\beta_{i} - {\frac{a_{i}\left( {\beta_{i + 1} - \beta_{i}} \right)}{a_{i + 1} - a_{i\;}}.}}}$

Then (24) becomes:

$\begin{matrix} {{\int\limits_{0}^{A}{{g(x)}{\mathbb{d}x}}} = {{{\sum\limits_{i = 1}^{N}{\int\limits_{a_{i}}^{a_{i + 1}}{m_{i}x}}} + {b_{i}{\mathbb{d}x}}} = {{{{\sum\limits_{i = 1}^{N}{m_{i}\frac{x^{2}}{2}}} + {b_{i}x}}❘_{a_{i}}^{a_{i + 1}}} = {{{\sum\limits_{i = 1}^{N}{m_{i}\left( {\frac{a_{i + 1}^{2}}{2} - \frac{a_{i}^{2}}{2}} \right)}} + {b_{i}\left( {a_{i + 1} - a_{i\;}} \right)}} = {{\sum\limits_{i = 1}^{N}{\left( \frac{\beta_{i + 1} - \beta_{i}}{a_{i + 1} - a_{i}} \right)\left( {\frac{a_{i + 1}^{2}}{2} - \frac{a_{i}^{2}}{2}} \right)}} + {\left\lbrack {\beta_{i} - \frac{a_{i}\left( {\beta_{i + 1} - \beta_{i}} \right)}{a_{i + 1} - a_{i}}} \right\rbrack\left( {a_{i + 1} - a_{i}} \right)}}}}}} & (25) \end{matrix}$

After further simplification, (25) becomes:

${\int\limits_{0}^{A}{{g(x)}{\mathbb{d}x}}} = {{\sum\limits_{i = 1}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}$

Where:

$\begin{matrix} {T_{i}^{1} = {a_{i + 1} - \frac{a_{i + 1}^{2} - a_{i}^{2}}{2\left( {a_{i + 1} - a_{i}} \right)}}} & (26) \\ {T_{i}^{2} = {\frac{a_{i + 1}^{2} - a_{i}^{2}}{2\left( {a_{i + 1} - a_{i}} \right)} - a_{i}}} & (27) \end{matrix}$

Define:

$\begin{matrix} {{g_{1}\left( {\beta_{1},{\ldots\mspace{14mu} B_{N + 1}}} \right)} = {{\sum\limits_{i = 1}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}} & (28) \end{matrix}$

Then the unity c.d.f. constraint, (4), is expressed as g ₁(β₁, . . . β_(N+1))=1  (29) 1.4. Constant Power Constraint

Consider (5):

${\int_{0}^{A}{x^{2}{g(x)}\ {\mathbb{d}x}}} = {{\int_{0}^{A}{x^{2}{\sum\limits_{i = 1}^{N}{{U_{i\;}(x)}\ {\mathbb{d}x}}}}} = {{\sum\limits_{i = 1}^{N}{\int_{0}^{A}{x^{2}{U_{i\;}(x)}\ {\mathbb{d}x}}}} = {{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{x^{2}{U_{i\;}(x)}\ {\mathbb{d}x}}}} = {{\sum\limits_{i = 1}^{N}{\int_{a_{i}}^{a_{i + 1}}{{x^{2}\left( {{m_{i}x} + b_{i}} \right)}{\mathbb{d}x}}}} = {{{{\sum\limits_{i = 1}^{N}{m_{i}\frac{x^{4}}{4}}} + {b_{i}\frac{x^{3}}{3}}}|_{a_{i}}^{a_{i + 1}}} = {{\sum\limits_{i = 1}^{N}{m_{i}\left( {\frac{a_{i + 1}^{4}}{4} - \frac{a_{i}^{4}}{4}} \right)}} + {b_{i}\left( {\frac{a_{i + 1}^{3}}{3} - \frac{a_{i}^{3}}{3}} \right)}}}}}}}$

This equation comes directly from equation 5, where we replaced the g(x) with the assumed form shown earlier. After much simplification, this equation becomes:

${\int_{0}^{A}{x^{2}{g(x)}{\mathbb{d}x}}} = {{\sum\limits_{i = 1}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}}}$

Where:

$W_{i}^{1} = {\frac{a_{i + 1}^{3} - a_{i}^{3}}{3} - \frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} + \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}}$ $W_{i}^{2} = {\frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} - \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}}$

Define:

$\begin{matrix} {{g_{2}\left( {B_{1},{\ldots\mspace{14mu}\beta_{N + 1}}} \right)} = {{\sum\limits_{i = 1}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}}}} & (30) \end{matrix}$

Then, the constant power constraint (5) becomes: g ₂(β₁, . . . β_(N+1))=σ². 1.5. Lagrangian Definition

When we take the partial derivatives of the Lagrangian and set them to zero to find the stationary points, we get a system of equations which includes the original system of equations (i.e., for F), plus the two constraint equations. Hence, the Lagrangian contains or includes or encompasses the constraints along with the original system, allowing the Lagrangian, i.e. the mathematical function that summarizes the dynamics of the system, to be defined as: Λ(β₁, . . . β_(N+1),λ₁,λ₂)=F(β₁, . . . β_(N+1))+λ₁ [g ₁(β₁, . . . β_(N+1))−1]+λ₂ [g ₂(β₁, . . . β_(N+1))−σ²]

Solve the set of equations:

$\begin{matrix} {{\frac{\partial\Lambda}{\partial\beta_{j}} = 0}{\frac{\partial\Lambda}{\partial\lambda_{1}} = 0}{\frac{\partial\Lambda}{\partial\lambda_{2}} = 0.}} & (31) \end{matrix}$

We have:

$\begin{matrix} {{\frac{\partial\Lambda}{\partial\beta_{j}} = {\frac{\partial F}{\partial\beta_{j}} + {\lambda_{1}\frac{\partial g_{1}}{\partial\beta_{j}}} + {\lambda_{2}\frac{\partial g_{2}}{\partial\beta_{j}}}}}{\frac{\partial\Lambda}{\partial\lambda_{1}} = {{g_{1}\left( {\beta_{1},{\ldots\mspace{14mu} B_{N + 1}}} \right)} - 1}}{\frac{\partial\Lambda}{\partial\lambda_{2}} = {{g_{2}\left( {\beta_{1},{\ldots\mspace{14mu} B_{N + 1}}} \right)} - {\sigma^{2}.}}}} & (32) \end{matrix}$

From (32), we see that we need the partial derivatives of g₁,g₂. We set the partial derivatives to zero, then solve for the betas to give values that under certain conditions represent the maximum or minimum of the function. We do this for the Lagrangian, because we want to include the constraints also, so we solve for beta parameter values which maximize or minimize the system and which also satisfy the constraints. From (28), we have, for j=1:

$\begin{matrix} {\frac{\partial g_{1}}{\partial\beta_{j}} = \left\{ \begin{matrix} T_{1}^{1} & {{{for}\mspace{14mu} j} = 1} \\ {T_{j - 1}^{2} + T_{j}^{1}} & {{{for}\mspace{14mu} 1} < j \leq N} \\ T_{N}^{2} & {{{for}\mspace{14mu} j} = {N + 1}} \end{matrix} \right.} & (33) \end{matrix}$

From (30), we have:

$\begin{matrix} {\frac{\partial g_{2}}{\partial\beta_{j}} = \left\{ \begin{matrix} W_{1}^{1} & {{{{for}\mspace{14mu} j} = 1}\;} \\ {W_{j - 1}^{2} + W_{j}^{1}} & {{{for}\mspace{14mu} 1} < j \leq N} \\ W_{N}^{2} & {{{for}\mspace{14mu} j} = {N + 1}} \end{matrix} \right.} & (34) \end{matrix}$

Combining all the partial derivative expressions from (23), (33) and (34) gives:

$\begin{matrix} {\frac{\partial\Lambda}{\partial\beta_{j}} = \left\{ \begin{matrix} {\left( {P_{0} + {P_{1}\beta_{1}} + {P_{2}\beta_{2}}} \right) + {\lambda_{1}T_{1}^{1}} + {\lambda_{2}W_{1}^{1}}} & {{{for}\mspace{14mu} j} = 1} \\ \begin{matrix} {\left\lbrack {{Q_{0}(j)} + {{Q_{1}(j)}\beta_{j - 1}} + {{Q_{2}(j)}\beta_{j}} + {{Q_{3}(j)}\beta_{j + 1}}} \right\rbrack +} \\ {{\lambda_{1}\left( {T_{j - 1}^{2} + T_{j}^{1}} \right)} + {\lambda_{2}\left( {W_{j - 1}^{2} + W_{j}^{1}} \right)}} \end{matrix} & {{{for}\mspace{14mu} 1} < j \leq N} \\ {\left( {R_{0\;} + {R_{1}\beta_{N}} + {R_{2}\beta_{N + 1}}} \right) + {\lambda_{1}T_{N}^{2}} + {\lambda_{2}W_{N}^{2}}} & {{{for}\mspace{14mu} j} = {N + 1}} \end{matrix} \right.} & (35) \end{matrix}$

From (35) and the system (31), we have the system of equations:

  (T₁¹)λ + (W₁¹)λ₂ + (P₁)β₁ + (P₂)β₂ = −P₀      ⋮(T_(j − 1)² + T_(j)¹)λ₁ + (W_(j − 1)² + W_(j)¹)λ₂ + (Q₁(j))β_(j − 1) + (Q₂(j))β_(j) + (Q₃(j))β_(j + 1) = −Q₀(j)      ⋮  (T_(N)²)λ₁ + (W_(N)²)λ₂ + (R₁)β_(N) + (R₂)β_(N + 1) = −R₀(0)λ₁ + (0)λ₂ + T₁¹β₁ + (T₁² + T₂¹)β₂ + … + (T_(N − 1)² + T_(N)¹)β_(N) + (T_(N)²)β_(N + 1) = 1 (0)λ₁ + (0)λ₂ + W₁¹β₁ + (W₁² + W₂¹)β₂ + … + (W_(N − 1)² + W_(N)¹)β_(N) + (W_(N)²)β_(N + 1) = σ².

The vertical three dots, as used above, signify that you continue on with the same pattern for the next value of the index j. So, the first equation gives the expression for j=1, then you continue with the same pattern for the remaining j values until the next expression. The reason why the three dots are used is two-fold: first it is compact notation, and second you don't have an exact number for the maximum value of j (it is a variable) so you can't write out all of the equations, in general.

In matrix form, the system of equations becomes:

$\begin{matrix} {{\begin{pmatrix} T_{1}^{1} & W_{1}^{1} & P_{1} & P_{2} & 0 & \; & \; & \; & \; & 0 \\ {T_{1}^{2} + T_{2}^{1}} & {W_{1}^{2} + W_{2}^{1}} & {Q_{1}(2)} & {Q_{2}(2)} & {Q_{3}(2)} & 0 & \; & \; & \; & 0 \\ \vdots & \; & \; & \; & \; & \; & \; & \; & \; & \; \\ {T_{j - 1}^{2} + T_{j}^{1}} & {W_{j - 1}^{2} + W_{j}^{1}} & 0 & \; & 0 & {Q_{1}(j)} & {Q_{2}(j)} & {Q_{3}(j)} & 0 & 0 \\ \vdots & \; & \; & \; & \; & \; & \; & \; & \; & \; \\ T_{N}^{2} & W_{N}^{2} & 0 & \; & \; & \; & \; & 0 & R_{1} & R_{2} \\ 0 & 0 & T_{1}^{1} & {T_{1}^{2} + T_{2}^{1}} & {T_{2}^{2} + T_{3}^{1}} & \ldots & \; & \; & {T_{N - 1}^{2} + T_{N}^{1}} & T_{N}^{2} \\ 0 & 0 & W_{1}^{1} & {W_{1}^{2} + W_{2}^{1}} & {W_{2}^{2} + W_{3}^{1}} & \ldots & \; & \; & {W_{N - 1}^{2} + W_{N}^{1}} & W_{N}^{2} \end{pmatrix}\begin{pmatrix} \lambda_{1} \\ \lambda_{2} \\ \; \\ \beta_{j} \\ \; \\ \beta_{N - 1} \\ \beta_{N} \\ \beta_{N + 1} \end{pmatrix}} = \begin{pmatrix} {- P_{0}} \\ {- {Q_{0}(2)}} \\ \vdots \\ {- {Q_{0}(j)}} \\ \vdots \\ {- R_{0}} \\ 1 \\ \sigma^{2} \end{pmatrix}} & (36) \end{matrix}$

The matrix form can then be solved numerically.

1.6. Hessian Matrix

A Hessian Matrix, sometimes referred to as simply a Hessian, is a square matrix of second-order partial derivatives of a scalar-valued function, or scalar field. It describes the local curvature of a function of many variables. This was also previously known by the term “functional determinants”.

In this case, a Hessian Matrix is useful in determination of the critical point type, which, after solving the system (36), can be obtained from the sign test on the sequence of principal minors of the bordered Hessian [23]. To calculate the Hessian, we need the second partial derivatives. For the Lagrangian system described above, the Hessian becomes the bordered Hessian, calculated as:

$H = \begin{pmatrix} 0 & 0 & \frac{\partial g_{1}}{\partial\beta_{1}} & \ldots & \frac{\partial g_{1}}{\partial\beta_{j}} & \ldots & \frac{\partial g_{1}}{\partial\beta_{N + 1}} \\ 0 & 0 & \frac{\partial g_{2}}{\partial\beta_{1}} & \ldots & \frac{\partial g_{2}}{\partial\beta_{j}} & \ldots & \frac{\partial g_{2}}{\partial\beta_{N + 1}} \\ \; & \; & \ddots & \; & \; & \; & \; \\ \; & \; & \; & \ddots & \; & {\frac{\partial^{2}F}{{\partial\beta_{j}}{\partial\beta_{i}}} + {\lambda_{1}\frac{\partial^{2}g_{1}}{{\partial\beta_{j}}{\partial\beta_{i}}}} + {\lambda_{2}\frac{\partial^{2}g_{2}}{{\partial\beta_{j}}{\partial\beta_{i}}}}} & \; \\ \; & \; & \; & \; & \ddots & \; & \; \end{pmatrix}$

From (33) and (34), we see that:

$\frac{\partial^{2}g_{1}}{{\partial\beta_{j}}{\partial\beta_{i}}} = {\frac{\partial^{2}g_{2}}{{\partial\beta_{j}}{\partial\beta_{i}}} = 0}$

From (23) we see that:

${\frac{\partial^{2}F}{\partial\beta_{1}^{2}} = P_{1}},{\frac{\partial^{2}F}{{\partial\beta_{2}}{\partial\beta_{1}}} = P_{2}},{\frac{\partial^{2}F}{{\partial\beta_{j}}{\partial\beta_{1}}} = {{0\mspace{14mu}{for}\mspace{14mu} j} > 2}}$ ${\frac{\partial^{2}F}{\partial\beta_{j}^{2}} = {Q_{2}(j)}},{\frac{\partial^{2}F}{{\partial\beta_{j - 1}}{\partial\beta_{j}}} = {Q_{1}(j)}},{\frac{\partial^{2}F}{{\partial\beta_{j + 1}}{\partial\beta_{j}}} = {Q_{3}(j)}},{\frac{\partial^{2}F}{{\partial\beta_{k}}{\partial\beta_{j}}} = {{0\mspace{14mu}{for}\mspace{14mu} k} \neq \left\{ {{j - 1},j,{j + 1}} \right\}}},{\frac{\partial^{2}F}{\partial\beta_{N + 1}^{2}} = R_{2}},{\frac{\partial^{2}F}{{\partial\beta_{N}}{\partial\beta_{N + 1}}} = R_{1}},{\frac{\partial^{2}F}{{\partial\beta_{j}}{\partial\beta_{N + 1}}} = {{0\mspace{14mu}{for}\mspace{14mu} j} < N}}$

Performing the sign test numerically shows that the solution to (36) is a minimum.

1.7 Numerical Solutions and Analysis

We now can solve system (36) numerically. Throughout, we take σ=√{square root over (2)}/2 to generate a normalized complex signal with unit power. FIGS. 2 through 9 contain plots of the optimal solution compared against the baseline Rayleigh distribution, as a function of cutoff amplitude value A. Values of N range from 8 to 11, depending on A. Above these N values, numerical problems occur causing the determinant of the system to be close to zero and thus precluding numerical solution.

Two observations are clear from the plots in FIGS. 2 through 9. The first observation is that as the value of A increases the solution approaches the baseline Rayleigh distribution. This is most apparent when comparing the solution prior to the tail of the Rayleigh distribution, since most of the solutions appear to match the Rayleigh distribution reasonably well before the tail. On the tail, the solution tends to curve upward to satisfy the unity c.d.f. and constant power constraints. This behavior suggests that closely matching the Rayleigh distribution when possible, is a desirable solution property.

The second observation is that for almost all of the solutions, the first ordinate, i.e. the y-coordinate, representing the distance from a point to the horizontal or x-axis measured parallel to the vertical or y axis, in this case β₁, is negative and thus a portion of the solution has negative values. This condition violates the non-negativity constraint on probability density functions.

There are two methods for overcoming the partial negativity of the optimal solution. The first method increases the first abscissa, the perpendicular distance of a point from the vertical axis, a₁=0 value to a small but positive value, specifically until β₁ in the solution becomes positive. In this case, the solution consists of the piecewise solution for x≧a₁, and the up-front portion of the Rayleigh distribution for 0≦x<a₁. The perturbation of a₁ has generally been observed to be small.

As an example, FIG. 10 contains the solution for A=1.4. In that instance, a value of a₁=0.0082 is sufficient to make β₁>0. The solution shown in FIG. 10 is negligibly different from the optimal solution shown in FIG. 5.

The second method for overcoming the negativity of the solution is to use the Rayleigh distribution as the front segment of the distribution, and then use a perturbed version of the back-end of the optimal solution. In this approach, the tail of the optimal solution must be perturbed to meet the new unity c.d.f. and constant power constraints, which are generated by replacing the front-end of the optimal solution with the Rayleigh distribution. This perturbation approach is discussed in the next section.

1.8. Perturbations from the Optimal Solution

Solutions of the system (36) tend to give a negative value for the first ordinate value β₁. Therefore, to determine a solution for which the non-negativity condition (6) holds, we perturb the optimal solution. We choose a_(k) starting at the distribution tail to perturb the solution. For xε[0,a₁] we use the original Rayleigh distribution ƒ_(R)(x), and we perturb the ordinate values for a_(j) with j>k for some k<N as depicted in FIG. 11.

We now define a perturbed function as:

$\begin{matrix} {{g_{ɛ}(x)} = {\sum\limits_{i = k}^{N}{U_{i}\left( {x,ɛ_{i}} \right)}}} & (37) \end{matrix}$

We define the perturbed function this way, because we intend to change the beta values by small amounts, i.e. by the epsilon values.

For some 1<k<N, where the U_(i)(x,ε_(i)) are perturbed, compactly-supported linear segments:

$\begin{matrix} {{U_{i}\left( {x,ɛ_{i}} \right)} = \left\{ \begin{matrix} {{\left( \frac{\left( {\beta_{i + 1} + ɛ_{i + 1}} \right) - \left( {\beta_{i} + ɛ_{i}} \right)}{a_{i + 1} - a_{i}} \right)\left( {x - a_{i}} \right)} + \left( {\beta_{i} + ɛ_{i}} \right)} & {{{for}\mspace{14mu} x} \in \left\lbrack {a_{i},a_{i + 1}} \right\rbrack} \\ 0 & {otherwise} \end{matrix} \right.} & (38) \end{matrix}$

Here, k represents the index value of the a value immediately before the a values corresponding to the beta values that we are going to change, see FIG. 11. We pick an index k, past which we adjust the beta parameters.

We then take ε_(k)=0 so that (a_(k), β_(k)+ε_(k))=(a_(k),β_(k)) thus leaving the point at index k fixed, maintaining a continuous function. Unity c.d.f and constant power constraints must still be satisfied for the perturbed solution; these are examined next.

1.8.1 Unity C.D.F. Constraint

From the unity c.d.f. condition on g_(ε)(x), we have:

1 = ∫₀^(A)g_(ɛ)(x) 𝕕x = ∫₀^(a_(k))f_(R)(x) 𝕕x + ∫_(a_(k))^(A)g_(ɛ)(x) 𝕕x

which gives:

$\begin{matrix} {{\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {1 - {\int_{0}^{a_{k}}{{f_{R}(x)}\ {\mathbb{d}x}}}}} & (39) \end{matrix}$

From integration of the Rayleigh distribution in (16), we obtain:

${\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {{1 - \left\lbrack {1 - {\mathbb{e}}^{- \frac{a_{k}^{2}}{\sigma^{2}}}} \right\rbrack} = {\mathbb{e}}^{- \frac{a_{k}^{2}}{\sigma^{2}}}}$

Now from (37):

${\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {{\int_{a_{k}}^{A}{\sum\limits_{i = k}^{N}{{U_{i}\left( {x,ɛ_{i}} \right)}\ {\mathbb{d}x}}}} = {{\sum\limits_{i = k}^{N}{\int_{a_{k}}^{A}{{U_{i}\left( {x,ɛ_{i}} \right)}\ {\mathbb{d}x}}}} = {{\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{{U_{i}\left( {x,ɛ_{i}} \right)}\ {\mathbb{d}x}}}} = {{\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{m_{i}x}}} + {b_{i}\ {\mathbb{d}x}}}}}}$

Where:

$m_{i} = \frac{\left( {\beta_{i + 1} - \beta_{i}} \right) + \left( {ɛ_{i + 1} - ɛ_{i}} \right)}{a_{i + 1} - a_{i}}$ $b_{i} = {\beta_{i} + ɛ_{i} - {\frac{a_{i}\left\lbrack {\left( {\beta_{i + 1} - \beta_{i}} \right) + \left( {ɛ_{i + 1} - ɛ_{i}} \right)} \right\rbrack}{a_{i + 1} - a_{i}}.}}$

Thus:

${{\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{m_{i}x}}} + {b_{i}\ {\mathbb{d}x}}} = {{{{\sum\limits_{i = k}^{N}{m_{i}\frac{x^{2}}{2}}} + {b_{i}x}}|_{a_{i}}^{a_{i + 1}}} = {{\sum\limits_{i = k}^{N}{m_{i}\left( {\frac{a_{i + 1}^{2}}{2} - \frac{a_{i}^{2}}{2}} \right)}} + {b_{i}\left( {a_{i + 1} - a_{i}} \right)}}}$

Which, after much simplification, gives:

$\begin{matrix} {{\int_{a_{k}}^{A}{{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {{\sum\limits_{i = k}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}} + {ɛ_{i}\left( \frac{a_{i + 1} - a_{i}}{2} \right)} + {ɛ_{i + 1}\left( \frac{a_{i + 1} - a_{i}}{2} \right)}}} & (40) \end{matrix}$

for T_(i) ¹ and T_(i) ² defined in (26) and (27).

Hence from (39) and (40) we get:

$\begin{matrix} {{{\sum\limits_{i = k}^{N}{ɛ_{i}\left( \frac{a_{i + 1} - a_{i}}{2} \right)}} + {ɛ_{i + 1}\left( \frac{a_{i + 1} - a_{i}}{2} \right)}} = {1 - {\int_{0}^{a_{k}}{{f_{R}(x)}\ {\mathbb{d}x}}} - {\sum\limits_{i = k}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}} & (41) \end{matrix}$ 1.8.2 Constant Power Constraint

In general, we want to keep the average power level the same as for the original, uncompanded OFDM signal, while simultaneously reducing the peak power, so as to lower the PAPR value of the companded OFDM signal relative to the nominal OFDM signal, so a constant power constraint is used.

From the constant power constraint on g_(ε)(x), we have:

σ² = ∫₀^(A)x²g_(ɛ)(x) 𝕕x = ∫₀^(a_(k))x²f_(R)(x) 𝕕x + ∫_(a_(k))^(A)x²g_(ɛ)(x) 𝕕x

which gives:

$\begin{matrix} {{\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {\sigma^{2} - {\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {\mathbb{d}x}}}}} & (42) \end{matrix}$

Using integration by parts, we get:

${\int_{c}^{d}{x^{2}{f_{R}(x)}\ {\mathbb{d}x}}} = {{\left( {c^{2} + \sigma^{2}} \right){\mathbb{e}}^{- \frac{c^{2}}{\sigma^{2}}}} - {\left( {d^{2} + \sigma^{2}} \right){\mathbb{e}}^{\frac{d^{2}}{\sigma^{2}}}}}$

and so:

${\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {\mathbb{d}x}}} = {\sigma^{2} - {\left( {a_{k}^{2} + \sigma^{2}} \right){\mathbb{e}}^{- \frac{a_{k}^{2}}{\sigma^{2}}}}}$

from which follows, using (42):

${\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {\left( {a_{k}^{2} + \sigma^{2}} \right){\mathbb{e}}^{- \frac{a_{k}^{2}}{\sigma^{2}}}}$

Now:

${\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {{\int_{a_{k}}^{A}{x^{2}{\sum\limits_{i = k}^{N}\ {{U_{i}\left( {x,ɛ_{i}} \right)}{\mathbb{d}x}}}}} = {{\sum\limits_{i = k}^{N}{\int_{a_{k}}^{A}{x^{2}\ {U_{i}\left( {x,ɛ_{i}} \right)}{\mathbb{d}x}}}} = {{\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{x^{2}\ {U_{i}\left( {x,ɛ_{i}} \right)}{\mathbb{d}x}}}} = {\sum\limits_{i = k}^{N}{\int_{a_{i}}^{a_{i + 1}}{{x^{2}\left( {{m_{i}x} + b_{i}} \right)}\ {\mathbb{d}x}}}}}}}$

Which, after much manipulation, becomes:

$\begin{matrix} {{\int_{a_{k}}^{A}{x^{2}{g_{ɛ}(x)}\ {\mathbb{d}x}}} = {{\sum\limits_{i = k}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}} + {ɛ_{i}\eta_{i}^{1}} + {ɛ_{i + 1}\eta_{i}^{2}}}} & (43) \end{matrix}$

Where:

$\eta_{i}^{1} = {{\frac{a_{i + 1}^{3} - a_{i}^{3}}{3} - \frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} + \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}} = W_{i}^{1}}$ $\eta_{i}^{2} = {{\frac{a_{i + 1}^{4} - a_{i}^{4}}{4\left( {a_{i + 1} - a_{i}} \right)} - \frac{a_{i}\left( {a_{i + 1}^{3} - a_{i}^{3}} \right)}{3\left( {a_{i + 1} - a_{i}} \right)}} = {W_{i}^{2}.}}$

From (43) using (42), we get:

$\begin{matrix} {{{\sum\limits_{i = k}^{N}{ɛ_{i}\eta_{i}^{1}}} + {ɛ_{i + 1}\eta_{i}^{2}}} = {\sigma^{2} - {\int_{0}^{a_{k}}{x^{2}{f_{R}(x)}\ {\mathbb{d}x}}} - {\sum\limits_{i = k}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}\beta_{i + 1}}}} & (44) \end{matrix}$

Combining (41) and (44), we get the under-constrained system of equations:

$\begin{matrix} {{{\sum\limits_{i = k}^{N}{\xi_{i}ɛ_{i}}} + {\xi_{i}ɛ_{i + 1}}} = C_{1}} & (45) \\ {{{\sum\limits_{i = k}^{N}{\eta_{i}^{1}ɛ_{i}}} + {\eta_{i}^{2}ɛ_{i + 1}}} = C_{2}} & (46) \end{matrix}$

Where:

$\xi_{i} = \frac{a_{i + 1} - a_{i}}{2}$ $C_{1} = {1 - {\int_{0}^{a_{k}}{f_{R}\;(x)\ {\mathbb{d}x}}} - {\sum\limits_{i = k}^{N}{T_{i}^{1}\beta_{i}}} + {T_{i}^{2}\beta_{i + 1}}}$ $C_{2} = {\sigma^{2} - {\int_{0}^{a_{k}}{x^{2}f_{R}\;(x)\ {\mathbb{d}x}}} - {\sum\limits_{i = k}^{N}{W_{i}^{1}\beta_{i}}} + {W_{i}^{2}{\beta_{i + 1}.}}}$

Next, the optimal values of ε_(i) are found.

1.8.3 Optimal Perturbation Parameters

Our goal here is to find the minimal perturbation g_(ε)(x) away from g(x). To do so, we'll choose parameters (ε_(k+1), . . . , ε_(N)) to minimize:

$\begin{matrix} {{F_{ɛ}\left( {{ɛ_{{k + 1},}\mspace{11mu}\ldots}\mspace{14mu},ɛ_{N + 1}} \right)} = {\sum\limits_{i = {k + 1}}^{N + 1}ɛ_{i}^{2}}} & (47) \end{matrix}$

First we rewrite (45) and (46) as:

$\begin{matrix} {C_{1} = {{{\sum\limits_{i = k}^{N}{\xi_{i}ɛ_{i}}} + {\xi_{i}ɛ_{i + 1}}} = {{{\xi_{k}ɛ_{k}} + \left( {{\sum\limits_{i = {k + 1}}^{N - 1}{\xi_{i - 1}ɛ_{i}}} + {\xi_{i}ɛ_{i}}} \right) + \left( {{\xi_{N - 1}ɛ_{N}} + {\xi_{N}ɛ_{N}}} \right) + {\xi_{N}ɛ_{N + 1}}} = {{{\xi_{k}ɛ_{k}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}} = {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}}}}}} & (48) \end{matrix}$

because, by definition, we took ε_(k)=0. Similarly from (46):

$\begin{matrix} {C_{2} = {{{\sum\limits_{i = k}^{N}{\eta_{i}^{1}ɛ_{i}}} + {\eta_{i}^{2}ɛ_{i + 1}}} = {{{\eta_{k}^{1}ɛ_{k}} + \left( {{\sum\limits_{i = {k + 1}}^{N - 1}{\eta_{i - 1}^{2}ɛ_{i}}} + {\eta_{i}^{1}ɛ_{i}}} \right) + \left( {{\eta_{N - 1}^{2}ɛ_{N}} + {\eta_{N}^{1}ɛ_{N}}} \right) + {\eta_{N}^{2}ɛ_{N + 1}}} = {{{\eta_{k}^{1}ɛ_{k}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}} = {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}}}}}} & (49) \end{matrix}$

Now we define functions:

$\begin{matrix} {{h_{1}\left( {ɛ_{k + 1},\ldots\mspace{14mu},ɛ_{N + 1}} \right)} = {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}}} & (50) \\ {{h_{2}\left( {ɛ_{k + 1},\ldots\mspace{14mu},ɛ_{N + 1}} \right)} = {\left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}}} & (51) \end{matrix}$

Then from (48) and (49) we get two constraint equations: h ₁(ε_(k+1), . . . ,ε_(N+1))=C ₁ h ₂(ε_(k+1), . . . ,ε_(N+1))=C ₂

Next, we define the Lagrangian: Λ_(h)(ε_(k+1), . . . ,ε_(N+1),λ₁,λ₂)=F _(ε)(ε_(k+1), . . . ,ε_(N+1))+λ₁ [h ₁(ε_(k+1), . . . ,ε_(N+1))−C ₁]+λ₂ [h ₂(ε_(k+1), . . . ε_(N+1))−C ₂]

Now, we solve the set of equations:

$\begin{matrix} {{\frac{\partial\Lambda_{h}}{\partial ɛ_{j}} = 0}{\frac{\partial\Lambda_{h}}{\partial\lambda_{1}} = 0}{\frac{\partial\Lambda_{h}}{\partial\lambda_{2}} = 0.}} & (52) \end{matrix}$

We have:

$\frac{\partial\Lambda_{h}}{\partial ɛ_{j}} = {\frac{\partial F_{ɛ}}{\partial ɛ_{j}} + {\lambda_{1}\frac{\partial h_{1}}{\partial ɛ_{j}}} + {\lambda_{2}\frac{\partial h_{2}}{\partial ɛ_{j}}}}$ $\frac{\partial\Lambda_{h}}{\partial\lambda_{1}} = {{h_{1}\left( {ɛ_{k + 1},\ldots\mspace{14mu},ɛ_{N + 1}} \right)} - C_{1}}$ $\frac{\partial\Lambda_{h}}{\partial\lambda_{2}} = {{h_{2}\left( {ɛ_{k + 1},\ldots\mspace{14mu},ɛ_{N + 1}} \right)} - C_{2}}$

From (47):

$\frac{\partial F_{ɛ}}{\partial ɛ_{j}} = {2ɛ_{j}}$

From (50) and (51):

$\frac{\partial h_{1}}{\partial ɛ_{i}} = \left\{ {\begin{matrix} {\xi_{i - 1} + \xi_{i}} & {{{{if}\mspace{14mu} k} + 1} \leq i \leq N} \\ \xi_{N} & {{{if}\mspace{14mu} i} = {N + 1}} \end{matrix},{\frac{\partial h_{2}}{\partial ɛ_{i}} = \left\{ \begin{matrix} {\eta_{i - 1}^{2} + \eta_{i}^{1}} & {{{{if}\mspace{14mu} k} + 1} \leq i \leq N} \\ \eta_{N}^{2} & {{{if}\mspace{14mu} i} = {N + 1}} \end{matrix} \right.}} \right.$

Thus:

$\frac{\partial\Lambda_{j}}{\partial ɛ_{j}} = \left\{ \begin{matrix} {{2ɛ_{i}} + {\lambda_{1}\left( {\xi_{i - 1} + \xi_{i\;}} \right)} + {\lambda_{2}\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)}} & {{{{if}\mspace{14mu} k} + 1} \leq i \leq N} \\ {{2ɛ_{N + 1}} + {\lambda_{1}\left( \xi_{N} \right)} + {\lambda_{2}\left( \eta_{N}^{2} \right)}} & {{{if}\mspace{14mu} i} = {N + 1}} \end{matrix} \right.$

Hence, we need to solve the following system of equations for ε_(i), λ₁, λ₂:

  (ξ_(k) + ξ_(k + 1))λ₁ + (η_(k)² + η_(k + 1)¹)λ₂ + (2)ɛ_(k + 1) = 0      ⋮  (ξ_(i − 1) + ξ_(i))λ₁ + (η_(i − 1)² + η_(i)¹)λ₂ + (2)ɛ_(i) = 0      ⋮  (ξ_(N))λ₁ + (η_(N)²)λ₂ + (2)ɛ_(N + 1) = 0 ${{(0)\lambda_{1}} + {(0)\lambda_{2}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\xi_{i - 1} + \xi_{i}} \right)ɛ_{i\;}}} \right) + {\left( {\xi_{N - 1} + \xi_{N}} \right)ɛ_{N}} + {\xi_{N}ɛ_{N + 1}}} = {{{{C_{1}(0)}\lambda_{1}} + {(0)\lambda_{2}} + \left( {\sum\limits_{i = {k + 1}}^{N - 1}{\left( {\eta_{i - 1}^{2} + \eta_{i}^{1}} \right)ɛ_{i\;}}} \right) + {\left( {\eta_{N - 1}^{2} + \eta_{N}^{1}} \right)ɛ_{N}} + {\eta_{N}^{2}ɛ_{N + 1}}} = C_{2}}$

Or:

$\begin{matrix} {{\begin{pmatrix} \left( {\xi_{k} + \xi_{k + 1}} \right) & \left( {\eta_{k}^{2} + \eta_{k + 1}^{1}} \right) & 2 & 0 & \; & \; & \; & \; & 0 \\ \left( {\xi_{k + 1} + \xi_{k + 2}} \right) & \left( {\eta_{k + 1}^{2} + \eta_{k + 2}^{1}} \right) & 0 & 2 & 0 & \; & \; & \; & 0 \\ \vdots & \vdots & \; & \; & \ddots & \; & \; & \; & \; \\ \; & \; & \; & \; & \; & \ddots & \; & \; & \; \\ \left( {\xi_{i + 1} + \xi_{i}} \right) & \left( {\eta_{i + 1}^{2} + \eta_{i}^{1}} \right) & 0 & \; & \; & \; & 2 & 0 & 0 \\ \vdots & \vdots & \; & \; & \; & \; & \; & \ddots & \; \\ \xi_{N} & \eta_{N}^{2} & 0 & \; & \; & \; & \; & 0 & 2 \\ 0 & 0 & \left( {\xi_{k} + \xi_{k + 1}} \right) & \; & \; & \ldots & \; & \; & \xi_{N} \\ 0 & 0 & \left( {\eta_{k}^{2} + \eta_{k + 1}^{1}} \right) & \; & \; & \ldots & \; & \; & \eta_{N}^{2} \end{pmatrix}\begin{pmatrix} \lambda_{1} \\ \lambda_{2} \\ \; \\ \; \\ {\; ɛ_{i}} \\ \; \\ ɛ_{N - 1} \\ ɛ_{N} \\ ɛ_{N + 1} \end{pmatrix}} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ \; \\ \; \\ \vdots \\ \; \\ C_{1} \\ C_{2} \end{pmatrix}} & (53) \end{matrix}$

This system of equations (53) may then be solved numerically.

1.9. Perturbed Solution Numerical Results

Perturbation solutions for the system in (53) are shown for A=1.4 and kε{5,6,7,8} in FIGS. 12 through 15. In these Figures, the perturbation solution with the Rayleigh front-end is plotted against the optimal solution.

1.10. Compander Derivation

Next, we denote the companding function by C, which operates on the input sequence x(n) by modifying the amplitude to produce the companded output sequence y(n)=C{x(n)}. The compander is monotonic, i.e. it maintains order of the function (i.e., either increasing or decreasing), over its domain xε[0,A], so the compander may be derived as: C{x(n)}=sgn{x(n)}F _(|y|) ⁻¹ F _(|x|) {x(n)}  (54)

Where F_(|x|) denotes the c.d.f of the input and F_(|y|) denotes the c.d.f of the output companded signal. To derive the compander, expressions for each c.d.f in (54) are needed. The p.d.f. of the companded signal is given by:

${f_{y}(x)} = \left\{ \begin{matrix} {\frac{2x}{\sigma^{2}}{\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}}} & {{{for}\mspace{14mu} 0} \leq x < a_{1}} \\ {f_{a_{1}}(x)} & {{{for}\mspace{14mu} a_{1}} \leq x < a_{2}} \\ \vdots & \vdots \\ {f_{a_{i}}(x)} & {{{for}\mspace{14mu} a_{i}} \leq x < a_{i + 1}} \\ \vdots & \vdots \\ {f_{a_{N}}(x)} & {{{for}\mspace{14mu} a_{N}} \leq x < a_{N + 1}} \\ 0 & {{{for}\mspace{14mu} a_{N + 1}} \leq x} \end{matrix} \right.$

For a_(i)≦x<a_(i+1), the c.d.f. is given by:

$\begin{matrix} {{F_{y}(x)} = {{\int_{0}^{x}{{f_{y}(\xi)}\ {\mathbb{d}\xi}}} = {\left\lbrack {{\int_{0}^{a_{1}}{{f_{R}(\xi)}\ {\mathbb{d}\xi}}} + {\sum\limits_{j = 1}^{i - 1}{\int_{a_{j}}^{a_{j + 1}}{{f_{a_{j}}\ (\xi)}{\mathbb{d}\xi}}}}} \right\rbrack + {\int_{a_{i}}^{x}{{f_{a_{i}}(\xi)}\ {\mathbb{d}\xi}}}}}} & (55) \\ {\mspace{79mu}{= {Z_{i - 1} + {\int_{a_{i}}^{x}{{f_{a_{i}}(\xi)}\ {\mathbb{d}\xi}}}}}} & (56) \end{matrix}$

Where Z_(i−1) denotes the term in brackets. Thus:

$\begin{matrix} {{{F_{y}(x)} - Z_{i - 1}} = {{\int_{0}^{x}{{f_{a_{i}}(\xi)}\ {\mathbb{d}\xi}}} = {{{\int_{a_{i}}^{x}{m_{i}\xi}} + {b_{i}\ {\mathbb{d}\xi}}} = {{{m_{i}\frac{x^{2}}{2}} + {b_{i}x} - {m_{i}\frac{a_{i}^{2}}{2}} - {b_{i}a_{i}}} = {{{\left( \frac{m_{i}}{2} \right)x^{2}} + {\left( b_{i} \right)x} + \left\lbrack {- \left( {{m_{i}\frac{a_{i}^{2}}{2}} + {b_{i}a_{i}}} \right)} \right\rbrack} = {{(\gamma)x^{2}} + {(\delta)x} + e}}}}}} & (57) \end{matrix}$

with the obvious definitions for γ, δ, e, i.e. they correspond to the coefficients of the x times x term, the x term, and the constant term. Now, γ might be negative, so divide both sides of (57) by γ to get:

$\begin{matrix} {\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} = {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \left( \frac{e}{\gamma} \right)}} & (58) \end{matrix}$

Complete the square in (58):

$\left( {x + \frac{\delta}{2\gamma}} \right)^{2} = {{x^{2} + {\frac{2\delta}{2\gamma}x} + \frac{\delta^{2}}{4\gamma^{2}}} = {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \left( \frac{\delta^{2}}{4\gamma^{2}} \right)}}$

Which implies:

$\begin{matrix} {{\left( {x + \frac{\delta}{2\gamma}} \right)^{2} + \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} = {{x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \frac{\delta^{2}}{4\gamma^{2}} + \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} = {x^{2} + {\left( \frac{\delta}{\gamma} \right)x} + \left( \frac{e}{\gamma} \right)}}} & (59) \end{matrix}$

Using (59), (58) becomes:

$\begin{matrix} {\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} = {\left. {\left( {x + \frac{\delta}{2\gamma}} \right)^{2} + \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}\Rightarrow\left( {x + \frac{\delta}{2\gamma}} \right)^{2} \right. = {\left. {\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}\Rightarrow{x + \frac{\delta}{\gamma^{2}}} \right. = \sqrt[ \pm ]{\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}}}} & (60) \end{matrix}$ (The above arrows are mathematical shorthand for the phrase: “which implies”)

And thus:

$\begin{matrix} {x = {\sqrt[ \pm ]{\frac{{F_{y}(x)} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}} & (61) \end{matrix}$

In (61) we need to decide whether to take the positive or negative square-root. The decision is made from (60) by examining the sign of

$x + {\frac{\delta}{2\gamma}.}$ If positive, the positive square-root is chosen, if negative, the negative square-root should be chosen. Of course, if x≧A, then we consider the sign on

$A + {\frac{\delta}{2\gamma}.}$ Finally, the compander is given by:

$\begin{matrix} {{F_{y}^{- 1}{F_{R}(x)}} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}} & (62) \end{matrix}$

for a_(i)≦x<a_(i+1). To complete the compander derivation, we need to calculate Z_(i) From (16) and (55):

$Z_{0} = {{\int_{0}^{a_{1}}{{f_{R}(\xi)}\ {\mathbb{d}\xi}}} = {1 - {\mathbb{e}}^{- \frac{a_{1}^{2}}{\sigma^{2}}}}}$

For i>0, Z_(i) contains integral terms of the form:

∫_(a_(j))^(a_(j + 1))f_(a_(j))(ξ) 𝕕ξ

which represents a trapezoidal area and can be shown to be:

${\int_{a_{j}}^{a_{j + 1}}{{f_{a_{j}}(\xi)}\ {\mathbb{d}\xi}}} = \frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}$

Hence:

$\begin{matrix} {Z_{i} = {Z_{i - 1} + \frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}}} & (63) \end{matrix}$ 1.11 Decompander Derivation

To find the decompander, we need to invert (62):

$\begin{matrix} {y = {\left. {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}\Rightarrow\left( {y + \frac{\delta}{2\gamma}} \right)^{2} \right. = {\left. {\frac{1 - {\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)}\Rightarrow{\gamma\left( {y + \frac{\delta}{2\gamma}} \right)}^{2} \right. = \left. {1 - {\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1} - e + \frac{\delta^{2}}{4\gamma}}\Rightarrow{{\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - {\gamma\left( {y + \frac{\delta}{2\gamma}} \right)}^{2} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}}\Rightarrow \right.}}} & (63) \\ {x = {{{sgn}(x)}\sqrt{{- \sigma^{2}}{\ln\left( {{- {\gamma\left( {y + \frac{\delta}{2\gamma}} \right)}^{2}} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}} \right)}}}} & (64) \end{matrix}$

Note that, in (64), the sgn function refers to retaining either the positive or negative square-root, depending on the sign of X.

Finally, we need to find the range of the compander for use in the decompander. When a_(i)≦x<a_(i+1), then y_(MIN)≦y<y_(MAX) where: y _(MIN)=min{y(a _(i)),y(a _(i+1))} y _(MAX)=max{y(a _(i)),y(a _(i+1))}

and from (63):

${y\left( a_{i} \right)} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{a_{i}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$ ${y\left( a_{i + 1} \right)} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{a_{i + 1}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$

2. SIMULATION

Simulation results were generated for a Quadrature Phase Shift Keying (QPSK) (a type of phase shift keying) modulated OFDM signal with sixty four subcarriers over an additive white Gaussian noise channel. Ten thousand realizations of signal-plus-noise were generated for each noise power.

For the Rayleigh p.d.f., σ=√{square root over (2)}/2 was chosen. An oversampling factor of four was used so that the discrete PAPR reasonably approximates the PAPR from a continuous system as is described Y. Wang, L.-H. Wang, J.-H. Ge, and B. Ai., “An Efficient Nonlinear Companding Transform for Reducing PAPR of OFDM Signals”, IEEE Trans. Broadcast., vol. 58, no. 4, pp.: 677-684, December 2012. PAPR reduction was measured using the complimentary cumulative distribution function (CCDF), which gives the probability that the PAPR is above a given threshold value. Symbol error rates, to measure demodulation performance, and power spectrums, to measure out-of-band power reduction, were also generated.

Lagrange compander results were then generated and compared against two state-of-the-art companders; the linear compander (LIN) in Y. Wang, L.-H. Wang, J.-H. Ge, and B. Ai., “An Efficient Nonlinear Companding Transform for Reducing PAPR of OFDM Signals”, IEEE Trans. Broadcast., vol. 58, no. 4, pp.: 677-684, December 2012, and the two-component, pieceweise linear compander (LIN2) in Y. Wang, J.-H. Ge, L.-H. Wang, J. Li, and B. Ai., “Nonlinear Companding Transform for Reduction of Peak-to-Average Power ratio in OFDM Systems”, IEEE Trans. Broadcast., vol. 59, no. 2, pp.: 369-375, June 2013. The LIN compander was state-of-the-art as of December 2012. The LIN compander was improved to create the LIN2 compander. A cutoff value of A=1.6281 was chosen, as this appears to be the maximum cutoff value obtainable with the LIN2 compander. The LIN compander can attain a maximum cutoff value of about A=1.44, and this value was used for the LIN compander results. Therefore the Lagrange compander solution can provide higher cutoff values than can be obtained with the LIN and LIN2 solutions, thereby providing additional compander design flexibility. Two Lagrange companders were used; the optimal solution (LG) with a perturbed a₁>0 to generate a non-negative solution, and the epsilon-perturbed compander (EP), which takes the optimal Lagrange compander and perturbs the tail using optimal perturbation values. For each plot, the curve labels have the following meaning: LG=optimal Lagrange, EP=epsilon-perturbed Lagrange, LW=LIN2, WG=LIN.

FIG. 16 contains symbol error plots for the four companders. FIG. 16 shows that the LG and EP companders provide an improvement in symbol error rate performance over the LIN and LIN2 companders.

This improvement in symbol error rate comes at the cost of a small reduction in PAPR performance relative to the LIN2 approach, shown in FIG. 17 which contains the CCDF results. The LIN approach, under the present scenario, gives the best performance due to the significantly reduced cutoff value, A=1.44 versus A=1.6281.

The Lagrange companders, however, provide a small performance improvement in out-of-band power rejection over the LIN2 and a significant performance advantage over the LIN approach (FIG. 18). The Lagrange approach shows an improvement of about 0.4 dB over the LIN2 approach and over 1.0 dB in out-of-band power rejection over the LIN approach.

Sample Lagrange companders, for different cutoff values, are shown in FIG. 19.

3. CONCLUSION

Those skilled in the art will appreciate that we have developed a family of companders for PAPR reduction in OFDM signals using a constrained optimization approach. We derived companders and decompanders and demonstrated that the new companders can provide performance improvements over current state-of-the art solutions. The new companders can also provide solutions over ranges of the cutoff values where the current state-of-the-art companders fail to exist. The set of companders developed in this paper increase the solution set of companders to tradeoff between demodulation performance, PAPR reduction, and out-of-band power rejection.

The foregoing description of the embodiments of the invention has been presented for the purposes of illustration and description. It is not intended to be exhaustive or to limit the invention to the precise form disclosed. Many modifications and variations are possible in light of this disclosure. It is intended that the scope of the invention be limited not by this detailed description, but rather by the claims appended hereto. 

What is claimed is:
 1. A signal transmitting apparatus for an orthogonal frequency division multiplexing (OFDM) system, comprising: a receiver, configured to receive and expand an orthogonal frequency-divided multiplexed signal in accordance with a dynamic range function; and a transmitter in communication with said receiver configured to compress the expanded orthogonal frequency-divided multiplexed signal and transmit the resulting compressed signal; wherein said receiver and said transmitter together comprise a compander and said dynamic range function comprises a conversion of a Rayleigh amplitude distribution, said Rayleigh amplitude distribution comprising a plurality of Rayleigh distributed signal amplitude values, to a piecewise linear distribution which optimally matches the Rayleigh amplitude distribution up to a cutoff amplitude value, and wherein said piecewise linear distribution is determined by defining and solving a constrained optimization problem using Lagrange multipliers.
 2. A method of companding an orthogonal frequency division multiplexed signal comprising: receiving and expanding an orthogonal frequency-divided multiplexed signal comprising a sufficiently large number of tones to form a Gaussian distributed signal, wherein components of said Gaussian distributed signal generate Rayleigh distributed signal amplitude values, wherein said Rayleigh distributed signal amplitude values comprise a Rayleigh amplitude distribution; converting said Rayleigh amplitude distribution to a piecewise linear distribution which optimally matches the Rayleigh amplitude distribution up to a cutoff amplitude value, and wherein said piecewise linear distribution is determined by defining and solving a constrained optimization problem using Lagrange multipliers; compressing the expanded orthogonal frequency-divided multiplexed signal using said constrained optimization problem; and transmitting the resulting compressed signal.
 3. The method of claim 2 wherein said Rayleigh amplitude distribution is decomposed into disjoint regions over which a plurality of amplitude weighting functions are used.
 4. The method of claim 2 wherein said constrained optimization approach comprises optimizing an objective function with respect to changeable elements in the presence of constraints on those elements.
 5. The method of claim 4 wherein said Rayleigh amplitude distribution is a probability distribution function and a first said constraint is that said Rayleigh amplitude distribution must integrate to unity.
 6. The method of claim 4 wherein a second said constraint is that power must be constant across the distribution.
 7. The method of claim 4 wherein an optimal solution is minimally perturbed to generate only non-negative solutions.
 8. The method of claim 4 wherein said Rayleigh amplitude distribution is defined as: ${f_{R\;}(x)} = {\frac{2x}{\sigma^{2}}{{\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}}.}}$
 9. The method of claim 8 wherein further comprising determining a function g(x) which minimizes: ∫₀^(A)(f_(R)(x) − g(x))² 𝕕x.
 10. The method of claim 9 wherein g(x) is chosen to be of piecewise linear form.
 11. The method of claim 10 wherein g(x) is defined as ${g(x)} = {\sum\limits_{i = 1}^{N}{{U_{i}(x)}.}}$
 12. The method of claim 2 wherein the dynamics of the system are defined as: Λ(β₁, . . . β_(N+1),λ₁,λ₂)=F(β₁, . . . β_(N+1))+λ₁ [g ₁(β₁, . . . β_(N+1))−1]+λ₂ [g ₂(β₁, . . . β_(N+1))−σ²].
 13. The method of claim 12 wherein a first constraint, that the cumulative distribution function must integrate to unity, is expressed as: g₁(β₁, . . . β_(N+1))=1.
 14. The method of claim 12 wherein a second constraint, that power must be constant, is expressed as: g₂(β₁, . . . β_(N+))=σ².
 15. The method of claim 2 wherein said constrained optimization problem can be defined as: Minimize: F(β₁, …  B_(N + 1)) = ∫₀^(A) (f_(R)(x) − g(x))²𝕕x for β_(i), subject to ∫₀^(A) g(x)𝕕x = 1  and ∫₀^(A) x²g(x)𝕕x = ∫₀^(∞)x²f_(R)(x) 𝕕x = σ², using Lagrange multipliers.
 16. The method of claim 2 wherein said compander is defined as: ${F_{y}^{- 1}{F_{R}(x)}} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{x^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$ for a_(i)≦x<a_(i+1), wherein $Z_{i} = {Z_{i - 1} + {\frac{{\beta_{i + 1}\left( {a_{i + 1} - a_{i}} \right)} + {\beta_{i}\left( {a_{i + 1} - a_{i}} \right)}}{2}.}}$
 17. The method of claim 16 wherein a decompander for use with said compander is defined as: ${x = {{{sgn}(x)}\sqrt{{- \sigma^{2}}{\ln\left( {{- {\gamma\left( {y + \frac{\delta}{2\gamma}} \right)}^{2}} - e + \frac{\delta^{2}}{4\gamma} + 1 - Z_{i - 1}} \right)}}}},$ where the sgn function refers to retaining either the positive or negative square-root, depending on the sign of x; and wherein the range of said compander can be defined as, when a_(i)≦x<a_(i+1), then y_(MIN)≦y<y_(MAX) where: y_(MIN)=min{y(a_(i)),y(a_(i+1))}, y_(MAX)=max{y(a_(i)), y(a_(i+1))}, subject to ${y\left( a_{i} \right)} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{a_{i}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - \frac{\delta}{2\gamma}}$ and ${y\left( a_{i + 1} \right)} = {\sqrt[ \pm ]{\frac{1 - {\mathbb{e}}^{- \frac{a_{i + 1}^{2}}{\sigma^{2}}} - Z_{i - 1}}{\gamma} - \left( {\frac{e}{\gamma} - \frac{\delta^{2}}{4\gamma^{2}}} \right)} - {\frac{\delta}{2\gamma}.}}$
 18. The method of claim 17 wherein a Hessian Matrix is used to determine critical point types. 